1.1.1 library Setup

library(rmarkdown)
library(rmdformats)
library(stargazer)
library(caret)
library(tidyverse)
library(DescTools)
library(car)
library(DT)
library(arm)
library(ggformula)
library(GGally)
library(ggm)
library(graphics)
library(Hmisc)
library(olsrr)
library(ggExtra)
library(PerformanceAnalytics)
library(readxl)
library(gvlma)
library(lawstat)

1.1.2 데이터 불러오고 정제하기

dong_num: categorical variable, “동”, “리” area: 아파트 면적 price: 아파트 가격(매매가, 단위: 천 원) floor: 아파트 층수 year: 아파트 건축 연도

층수, 건축 연도, 집 넓이와 마지막으로 아파트가 시가지에 위치해 있는지 여부는 모두 아파트의 가격과 연관성이 있는지 확인해볼 가치가 있는 변수다. 이에 원 데이터에서 위의 다섯개 변수를 가져왔고, 상관성 분석을 통해 변수 간 상관성을 확인해 다중회귀분석을 시행할 것이다.

원데이터에 있으나 분석에 사용하지 않은 변수는 주소(번지, 동명, 아파트명, 도로명주소)과 계약일이다.

결측치가 있는 행은 제거했으며, 아파트 거래가 아닌 것으로 확인된 두 개의 influential point 역시 잘못된 데이터이므로 배제하고 분석을 진행한다.(상세한 내용 참고: 뒷부분 influential point, influential.test 부분에서 결정)

# 불러오기
ht <- read_excel("data/housing transaction.xlsx", skip = 16)
ht <- ht[, c(1,6, 9, 10, 11)] # 뒤의 가정검토 부분의 outlier test 결과 참조, 아파트가 아닌 자료 포함돼 있음. 

colnames(ht) <- c("dong", "area", "price", "floor", "year")

# 동/리 분류 및 더미변수 생성
for (i in 1:320) {
  ht[i,1] <- str_replace(ht[i, 1], "^경기도 파주시 ..동$", "동")
  ht[i,1] <- str_replace(ht[i, 1], "^경기도 파주시 .*읍 .*리$", "리")
  # character인 집값 데이터 정리
  ht[i,3] <- str_replace(ht[i, 3], ",", "")
}

# 결측치 제거
ht <- na.omit(ht)

# 데이터 as. numeric
ht$price <- as.numeric(ht$price)
ht$dong <- as.factor(ht$dong)
ht$dong_num<- ifelse(ht$dong == "동", 1, 0)
ht$area<- as.numeric(ht$area)
ht$floor<- as.numeric(ht$floor)
ht$year <- as.numeric(ht$year)
ht$logprice <- log2(ht$price)

ht_1 <- ht
ht <- ht[-c(11,19),]

1.1.3 데이터 확인하기

str(ht)

## tibble [318 x 7] (S3: tbl_df/tbl/data.frame)
##  $ dong    : Factor w/ 2 levels "동","리": 1 1 1 1 1 1 1 1 1 1 ...
##  $ area    : num [1:318] 84.8 84.8 75 59.9 85 ...
##  $ price   : num [1:318] 14000 14800 14000 15000 18750 ...
##  $ floor   : num [1:318] 1 3 12 2 7 17 17 18 18 2 ...
##  $ year    : num [1:318] 2000 2000 2001 2000 2000 ...
##  $ dong_num: num [1:318] 1 1 1 1 1 1 1 1 1 1 ...
##  $ logprice: num [1:318] 13.8 13.9 13.8 13.9 14.2 ...

datatable(ht)

attach(ht)
# 매매가 분포
hist(price, breaks="FD", col="gray", freq=FALSE, ylab="Density")
lines(density(price))

종속변수로 다룰 매매가다. 매매가는 정규분포와 비교해 좌측으로 치우쳐진 분포를 보여준다. 모든 변수의 분포를 확인하면 다음과 같다.

# 시각적으로 전체 분포 확인하기
ggpairs(data = ht, columns = c("price", "year", "area", "floor", "dong"))

또한 더미변수인 동/리에 따라 year, floor, area 별 price분포도 달라진다는 것을 확인할 수 있다. 특히, 셋 모두 동 더미변수가 일정할 때의 price를 종속변수로 두는 추세선들이 평행하지 않아 dong 더미변수와의 interaction이 유의한지 확인해볼 필요가 있으며, 특히 year변수 X축에 둘 때 추세선이 inordinal하게 엇갈려 interaction effect이 있는지, 즉 아파트가 동에 위치하는지 리에 위치하는지에 따라 매매가와 아파트 연식 간의 연관성이 다르게 나타나지 않는지 주의를 기울일 필요가 있다.

attach(ht)
car::scatterplot(price ~ year | dong_num, data=ht, 
            boxplots=FALSE,
            span=0.75, 
            col=c("mistyrose1", "mistyrose4"),
            legend.coords = "left",
            main="Price by Year")

car::scatterplot(price ~ floor | dong_num, data=ht, 
            boxplots=FALSE,
            span=0.75, 
            col=c("mistyrose1", "mistyrose4"), 
            legend.coords = "left",
            main="Price by Floor")

car::scatterplot(price ~ area | dong_num, data=ht, 
            boxplots=FALSE,
            span=0.75, 
            col=c("mistyrose1", "mistyrose4"), 
            legend.coords = "left",
            main="Price by Area")

ggplot() + geom_boxplot(data = ht, aes(x = dong, y = price)) + theme_light() + labs(title = "Price by Dong")

dong_num을 각각 0과 1로 일정하게 둘 때, 각각의 x변수 year, floor, area에 대해 추세선이 다르게 나타나는 것을 scatterplot에서 확인할 수 있다. area의 경우 특정 값에 샘플이 몰려있는 점을 유의해야 할 것으로 보인다.

세 변수를 모두 predictor로 하는 회귀분석을 진행할 예정이다.

1.1.4 소결

year 건축 연도가 최신일수록, floor 층수가 높을수록, area 면적이 넓을수록, dong 리보다는 동일 때 price 더 높은 가격에 매매가 이뤄진다.

floor 변수의 경우 매우 고르게 분포해있는데 반해, area의 경우 특정 면적의 표본 수가 컸는데, 이는 아파트 면적과 형태가 어느정도 정형화돼 있기 때문일 것이다.

dong 동에서의 거래량이 리에서의 거래 건수보다 훨씬 많고, 특히 동에서 리에서보다 아파트 건축 연수에 따른 아파트 매매가 차이가 컸다. 이는 아파트 자체가 동에 더 많고, 거래량도 더 많으며, 택지 개발, 신도시 개발로 동 전환을 한 지역이 도시 내에 분포해있기 때문 등에서 이유를 찾을 수 있다.

1.2.1 Correlation Matrix

# marginal correlation matrix
(cor.m <- cor(data.frame(price, year, floor, area, dong_num)))

##              price      year      floor        area    dong_num
## price    1.0000000 0.7793920 0.27268511  0.43464821  0.31383179
## year     0.7793920 1.0000000 0.18101250  0.10727676  0.15915674
## floor    0.2726851 0.1810125 1.00000000  0.01645352  0.09266305
## area     0.4346482 0.1072768 0.01645352  1.00000000 -0.06924280
## dong_num 0.3138318 0.1591567 0.09266305 -0.06924280  1.00000000

(cor.log <- cor(data.frame(log2(price), year, floor, area, dong_num)))

##             log2.price.      year      floor        area    dong_num
## log2.price.   1.0000000 0.7958411 0.23125275  0.47601586  0.34405680
## year          0.7958411 1.0000000 0.18101250  0.10727676  0.15915674
## floor         0.2312527 0.1810125 1.00000000  0.01645352  0.09266305
## area          0.4760159 0.1072768 0.01645352  1.00000000 -0.06924280
## dong_num      0.3440568 0.1591567 0.09266305 -0.06924280  1.00000000

# partial correlation between the pairs given the rest
parcor(cor.m)

##              price       year       floor       area    dong_num
## price    1.0000000  0.8254002  0.25810025  0.6356138  0.43176983
## year     0.8254002  1.0000000 -0.12160655 -0.4732295 -0.27786336
## floor    0.2581002 -0.1216066  1.00000000 -0.1619878 -0.05414133
## area     0.6356138 -0.4732295 -0.16198779  1.0000000 -0.33566780
## dong_num 0.4317698 -0.2778634 -0.05414133 -0.3356678  1.00000000

stargazer(rcorr(as.matrix(data.frame(price, year, floor, area, dong_num)), type = "pearson")[1], rcorr(as.matrix(data.frame(price, year, floor, area, dong_num)), type = "pearson")[3], column.labels = c("Correlation b/w variables", "Correlation p-value"), type = 'html')

price변수와 year floor area dong_num은 모두 상관계수가 있고, 모두 유의해 다중회귀분석을 진행해보기로 한다.

partial correlation matrix의 dong_num과 유의수준 0.05하에서 유의한 year는 .16의 상관계수가 있는데, 앞서 언급한 대로 파주시의 2000년대 시가지 개발 및 택지 개발, 2010년대초 신도시 개발과 연관있는 것으로 판단된다.

때문에 interaction 항에도 주목할 필요가 있다.

1.2.2 ANOVA

RQ: 주택 가격은 동/리 여부에 따라 다를까?

Is price of an appartment differ across location(city or not)?

H0: 동/리 여부와 주택 가격은 상관 없다. H1: 동/리 여부와 주택 가격은 상관 없지 않다.

total.mean <- mean(price)
total.var <- var(price)
gmean <- aggregate(price~dong, data = ht, mean)
gvar <- aggregate(price~dong, data = ht, var)
gn <- aggregate(price~dong, data = ht, length)
data.frame(total.mean, total.var) 

##   total.mean total.var
## 1   27886.16 117725817

cbind(`표본평균` = gmean, `표본분산` = gvar, `표본수` = gn)

##   표본평균.dong 표본평균.price 표본분산.dong 표본분산.price 표본수.dong
## 1            동       29337.17            동      117738053          동
## 2            리       19920.41            리       43535721          리
##   표본수.price
## 1          269
## 2           49

mod.anova <- aov(price~dong, data = ht)
summary(mod.anova)

##              Df    Sum Sq   Mean Sq F value   Pr(>F)    
## dong          1 3.676e+09 3.676e+09   34.52 1.07e-08 ***
## Residuals   316 3.364e+10 1.065e+08                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.3.1 RQ

Research Question(Question of Interest): 아파트 건축 연도, 아파트의 넓이와 층수, 도시에 위치하는지 여부는 아파트 가격 변화를 설명할 수 있는가?

Is price of house associated with year built, area(m2), floor, location(city or not) of housing?

\(\hat{Price} = \beta_0 + \beta_1Year + \beta_2Area + \beta_3Floor + \beta_4dong\)

설명변수: year area floor, dong_num 종속변수: price

1.3.2 변수별 단순선형회귀분석

예측변수: dong 종속변수: logprice

mod.dong <- lm(logprice ~ dong, data = ht); summary(mod.dong)

## 
## Call:
## lm(formula = logprice ~ dong, data = ht)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.89572 -0.31591  0.03201  0.35394  1.19044 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.74113    0.03361 438.617  < 2e-16 ***
## dong리      -0.55769    0.08562  -6.514 2.89e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5512 on 316 degrees of freedom
## Multiple R-squared:  0.1184, Adjusted R-squared:  0.1156 
## F-statistic: 42.43 on 1 and 316 DF,  p-value: 2.886e-10

0.05 유의수준에서 dong은 logPrice를 유의하게 설명한다. 다만 logPrice에 대해 11.56%의 설명력을 가져 설명력은 그리 높지 않다.

예측변수: year area floor 종속변수: logprice

mod.year <- lm(logprice ~ year, data=ht)
mod.area <- lm(logprice ~ area, data=ht)
mod.floor <- lm(logprice ~ floor, data=ht)
summary(mod.year);

## 
## Call:
## lm(formula = logprice ~ year, data = ht)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.25117 -0.23838  0.04659  0.24042  1.05042 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.457e+02  6.863e+00  -21.23   <2e-16 ***
## year         7.993e-02  3.421e-03   23.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3555 on 316 degrees of freedom
## Multiple R-squared:  0.6334, Adjusted R-squared:  0.6322 
## F-statistic: 545.9 on 1 and 316 DF,  p-value: < 2.2e-16

새 아파트: 0.05 유의수준에서 year은 logPrice를 유의하게 설명한다. 또한 logPrice에 대해 63.2%의 설명력을 가진다.

summary(mod.area);

## 
## Call:
## lm(formula = logprice ~ area, data = ht)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0452 -0.2791 -0.0226  0.3531  1.2813 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.625855   0.110827 122.947   <2e-16 ***
## area         0.011865   0.001233   9.622   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5163 on 316 degrees of freedom
## Multiple R-squared:  0.2266, Adjusted R-squared:  0.2241 
## F-statistic: 92.58 on 1 and 316 DF,  p-value: < 2.2e-16

넓은 아파트: 0.05 유의수준에서 area은 logPrice를 유의하게 설명한다. 또한 logPrice에 대해 22.41%의 설명력을 가진다.

summary(mod.floor)

## 
## Call:
## lm(formula = logprice ~ floor, data = ht)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.25015 -0.34208  0.05709  0.42448  1.12483 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.426228   0.062946 229.186  < 2e-16 ***
## floor        0.022328   0.005284   4.225 3.12e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5711 on 316 degrees of freedom
## Multiple R-squared:  0.05348,    Adjusted R-squared:  0.05048 
## F-statistic: 17.85 on 1 and 316 DF,  p-value: 3.125e-05

높은 층의 아파트: 0.05 유의수준에서 floor은 logPrice를 유의하게 설명한다. 또한 logPrice에 대해 5.05%의 설명력을 가져 이번 분석에서 다룰 설명변수 중 가장 설명력이 낮다. 그러나 coefficient가 유의하기 때문에 배제하지 않고 다중선형회귀분석을 진행한다.

stargazer(mod.year, mod.dong, mod.area, mod.floor, title= "각 설명변수 별 단순회귀분석 결과 비교", align = T, type = "html", column.labels = c("새 아파트인가", "동에 있나", "넓은 집인가", "높은 층인가"))

1.4.1 RQ

Research Question(Question of Interest): 동/리 여부에 따른 아파트 건축 연도, 아파트의 넓이와 층수 아파트 가격 변화를 설명할 수 있는가?

예측변수: year area dong_num floor 종속변수: logprice

H0: \(\beta_1=\beta_2=\beta_3=\beta_4=0\) H1: \(\beta_1,\beta_2,\beta_3,\beta_4\) is not all \(0\)

회귀 모델:

\(log\hat{Price} = \beta_0 + \beta_1Year + \beta_2Area + \beta_3Floor + \beta_4dong\)

동일 때 \(dong = 1\) \(log\hat{Price} = (\beta_0 + \beta_4) + \beta_1Year + \beta_2Area + \beta_3Floor\)

리일 때 \(dong = 0\) \(log\hat{Price} = \beta_0 + \beta_1Year + \beta_2Area + \beta_3Floor\)

설명변수: year area floor, dong_num 종속변수: price

1.4.2 분석

lm.me <- lm(logprice ~ year + area + floor + dong, data = ht)
summary(lm.me)

## 
## Call:
## lm(formula = logprice ~ year + area + floor + dong, data = ht)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.69759 -0.12131  0.00674  0.13050  0.53685 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.267e+02  4.426e+00 -28.623  < 2e-16 ***
## year         6.999e-02  2.210e-03  31.671  < 2e-16 ***
## area         1.041e-02  5.351e-04  19.456  < 2e-16 ***
## floor        7.201e-03  2.092e-03   3.443 0.000655 ***
## dong리      -4.136e-01  3.512e-02 -11.776  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2219 on 313 degrees of freedom
## Multiple R-squared:  0.8585, Adjusted R-squared:  0.8567 
## F-statistic: 474.8 on 4 and 313 DF,  p-value: < 2.2e-16

Dong의 동, 리 여부에 따라서도 Coefficient를 계산한다.

coef.me <- data.frame(coef(lm.me))
b0 <- coef.me[1,1]
b1 <- coef.me[2,1]
b2 <- coef.me[3,1]
b3 <- coef.me[4,1]
b4 <- coef.me[5,1] 

#전체
data.frame(b0, b1, b2, b3, b4)

##          b0         b1         b2         b3         b4
## 1 -126.6736 0.06999171 0.01041114 0.00720134 -0.4135748

#동
data.frame(b0 + b4, b1, b2, b3)

##     b0...b4         b1         b2         b3
## 1 -127.0872 0.06999171 0.01041114 0.00720134

#리
data.frame(b0, b1, b2, b3)

##          b0         b1         b2         b3
## 1 -126.6736 0.06999171 0.01041114 0.00720134

1.4.3 회귀 분석 결과

\(log\hat{Price} = -126.674+ 0.070Year + 0.010Area + 0.007Floor - 0.414dong\)

모든 설명변수는 매우 유의한 p-value 값을 가지며,(유의수준 \(\alpha = 0.05\)) H0은 기각된다. 즉 모든 beta는 0이 아니다.

위의 다중회귀모형은 Adjusted \(R^2\) 역시 .857로, 약 85.7%의 설명력을 갖는다.

더 상세하게 보면, 다른 변수가 통제된 상태에서, Year 한 단위 변화는(1년 더 최근에 지어졌을수록) Price 약 7.25%의 변화에 정방향으로 영향을 받으며, Area \(1m^2\) 증가할 때마다, 또 Floor 한 층 높아질 때마다 각각 Price는 약 1.01%, 0.70%의 변화에 정방향으로 영향이 있고, dong 한 단위 변화(리)는 Price의 -51.29% 변화에 영향이 있다.

categorical 변수에 따른 sector별 분석도 가능하다.

동일 때 \(Dong = 1\) \(log\hat{Price} = -127.087 + 0.070Year + 0.010Area + 0.007Floor\)

리일 때 \(dong = 0\) \(log\hat{Price} = -126.673 + 0.067Year + 0.012Area + 0.008Floor\)

두 섹터 모두 Main effect model의 가정에 따라 다른 변수가 통제된 상태에서 Year Area Floor 각각 한 단위 증가할 때, Price의 변화율은 서로 동일하고, 섹터를 나누지 않을 때와도 일정하다.

그러나 \(\beta0\)에서 차이점을 보인다.

결론의 영향이 있다, 영향을 받다는 association을 설명하는 것일 뿐, 인과관계의 의미를 내포하지 않는다.

1.4.4 시각화

par(mfrow = c(1,3))
betas.year <- coef(lm(logprice~year+dong, data = ht))

intercept.bc <- betas.year[1]
slope.year <- betas.year[2]
intercept.dong <- betas.year[1] + betas.year[3]
intercept.li   <- betas.year[1]

plot(logprice ~ year, data=ht,
            pch=20,
            xlab="year",
            ylab="log2Price", 
            main="House Price by year\n(Main-Effect Model)")

curve(intercept.dong + slope.year*x, col="aquamarine3", add=TRUE)
curve(intercept.li+ slope.year*x, col="aquamarine4", add=TRUE)

legend("topleft", c("동", "리"), 
                  lty=1, lwd=2, 
                  col=c("aquamarine3", "aquamarine4"),
                  text.col=c("aquamarine3", "aquamarine4"))

betas.area <- coef(lm(logprice~area+dong, data = ht))

intercept.bc <- betas.area[1]
slope.area <- betas.area[2]
intercept.dong <- betas.area[1] + betas.area[3]
intercept.li   <- betas.area[1]

plot(logprice ~ area, data=ht,
            pch=20,
            xlab="area",
            ylab="log2Price", 
            main="House Price by area\n(Main-Effect Model)")

curve(intercept.dong + slope.area*x, col="aquamarine3", add=TRUE)
curve(intercept.li+ slope.area*x, col="aquamarine4", add=TRUE)

legend("topleft", c("동", "리"), 
                  lty=1, lwd=2, 
                  col=c("aquamarine3", "aquamarine4"),
                  text.col=c("aquamarine3", "aquamarine4"))


betas.floor <- coef(lm(logprice~floor+dong, data = ht))

intercept.bc <- betas.floor[1]
slope.floor <- betas.floor[2]
intercept.dong <- betas.floor[1] + betas.floor[3]
intercept.li   <- betas.floor[1]

plot(logprice ~ floor, data=ht,
            pch=20,
            xlab="floor",
            ylab="log2Price", 
            main="House Price by floor\n(Main-Effect Model)")

curve(intercept.dong + slope.floor*x, col="aquamarine3", add=TRUE)
curve(intercept.li+ slope.floor*x, col="aquamarine4", add=TRUE)

legend("topleft", c("동", "리"), 
                  lty=1, lwd=2, 
                  col=c("aquamarine3", "aquamarine4"),
                  text.col=c("aquamarine3", "aquamarine4"))

각각 동, 리일 때 year, area, floor의 회귀 계수 값은 일정하다. 따라서 동, 리 직선이 평행하게 나타난다.

1.5.1 Estimation formula

\(log\hat{Price} = \beta_0 + \beta_1Year + \beta_2Area + \beta_3Floor + \beta_4Dong + \beta_5Dong^*Year\)

동일 때 \(Dong = 1\) \(log\hat{Price} = \beta_0+\beta_4 + \beta_1+\beta_5Year + \beta_2Area + \beta_3Floor\)

리일 때 \(dong = 0\) \(log\hat{Price} = \beta_0 + \beta_1Year + \beta_2Area + \beta_3Floor\)

설명변수: year area floor, dong 종속변수: price

1.5.2 Scope of Inference

이 회귀 분석의 결과를 더 general 한 범위의 추론에 활용할 수 없다. 변인을 통제한 experiment가 아니라 observation에 대한 사후 분석자료고, 한정된 지역과 시기의 자료기 때문이다.

1.5.3 분석

lm.iem <- lm(logprice ~ year + area + floor + dong + dong*year, data = ht)
summary(lm.iem)

## 
## Call:
## lm(formula = logprice ~ year + area + floor + dong + dong * year, 
##     data = ht)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.84173 -0.12527 -0.00214  0.14334  0.53148 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.350e+02  4.859e+00 -27.791  < 2e-16 ***
## year         7.416e-02  2.425e-03  30.576  < 2e-16 ***
## area         1.045e-02  5.241e-04  19.937  < 2e-16 ***
## floor        6.701e-03  2.053e-03   3.265 0.001217 ** 
## dong리       3.969e+01  1.055e+01   3.763 0.000200 ***
## year:dong리 -2.000e-02  5.261e-03  -3.803 0.000172 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2173 on 312 degrees of freedom
## Multiple R-squared:  0.8648, Adjusted R-squared:  0.8626 
## F-statistic:   399 on 5 and 312 DF,  p-value: < 2.2e-16

year*dong_num 항(교차항)은 0.05 수준에서 유의한 통계적 상관관계가 존재한다. coefficient \(\beta_0= -135.0\), \(\beta_5= 0.020\) Interaction model은 \(R^2 = 0.8626\)로 86.3%의 설명력을 가져 Main effect model에서보다 조금 더 높은 설명력을 갖는다.

동과 리의 model을 모양대로 따로 회귀 분석을 시행해 나머지 \(\beta\) 값을 구할 수도 있다.

coef.iem <- data.frame(coef(lm.iem))
b0 <- coef.iem[1,1]
b1 <- coef.iem[2,1]
b2 <- coef.iem[3,1]
b3 <- coef.iem[4,1]
b4 <- coef.iem[5,1] 
b5 <- coef.iem[6,1]

data.frame(b0, b1, b2, b3, b4, b5)

##          b0         b1         b2          b3      b4        b5
## 1 -135.0292 0.07415685 0.01044834 0.006700931 39.6862 -0.020005

#dong
data.frame(b0 + b4, b1 + b5, b2, b3)

##     b0...b4    b1...b5         b2          b3
## 1 -95.34301 0.05415185 0.01044834 0.006700931

#li
data.frame(b0, b1, b2, b3)

##          b0         b1         b2          b3
## 1 -135.0292 0.07415685 0.01044834 0.006700931

1.5.4 회귀 분석 결과

\(log\hat{Price} = -135.029 + 0.074Year + 0.010Area + 0.007Floor + 39.686Dong - 0.020Dong^*Year\)

모든 설명변수는 매우 유의한 p-value 값을 가지며,(유의수준 \(\alpha = 0.05\)) H0은 기각된다. 즉 모든 beta는 0이 아니다.

위의 다중회귀모형은 Adjusted \(R^2\) 역시 .863로, 약 86.3%의 설명력을 갖는다.

더 상세하게 보면, 다른 변수가 통제된 상태에서, Area \(1m^2\) 증가할 때마다, 또 Floor 한 층 높아질 때마다 각각 Price는 약 1.01%, 0.70%의 변화에 정방향으로 영향이 있다.

Year은 Dong 변수와 유의수준 0.05 에서 유의하게 interact하며, 따라서 Year 변수의 계수는 아파트가 동에 위치했느냐, 리에 위치했느냐에 따라 달라진다.

categorical 변수에 따른 sector별 분석도 가능하다. Dong에 따라 각 변수의 coefficient도 다르다.

동일 때 \(Dong = 1\) \(log\hat{Price} = -93.343 + 0.054Year + 0.010Area + 0.007Floor\)

더 상세하게 보면, 다른 변수가 통제된 상태에서, 아파트가 1년 더 최근에 지어졌을수록, Price 약 5.55%의 변화에 정방향으로 영향을 받으며, Area \(1m^2\) 넓어질 때마다, 또 Floor 한 층 높아질 때마다 각각 Price는 약 1.01%, 0.70%의 변화에 정방향으로 영향이 있다.

리일 때 \(dong = 0\) \(log\hat{Price} = -135.029 + 0.007Year + 0.010Area + 0.006Floor\)

\(\beta_0\)는 아파트가 동에 있을 때보다 작아 아파트의 매매가가 리에서 전반적으로 낮음을 알 수 있다.

더 상세하게 보면, 다른 변수가 통제된 상태에서, 아파트가 1년 더 최근에 지어졌을수록, Price 약 0.70%의 변화에 정방향으로 영향을 받으며, Area \(1m^2\) 넓어질 때마다, 또 Floor 한 층 높아질 때마다 각각 Price는 약 1.01%, 0.60%의 변화에 정방향으로 영향이 있다.

두 섹터에서 주목할 부분은 Year의 회귀 계수가 매우 다르다는 점이다. 특히 아파트가 최근에 지어졌는지 여부는 아파트가 동에 있을 때 계수 값이 커지는 것으로 드러났다. 이는 동지역일 경우 최근에 가까워질수록 역 인근 택지나 대형 아파트 단지가 더 지어졌기 때문일 것으로 추정된다.

\(\beta0\)에서도 차이점이 보여, 동지역이 전반적으로 매매가가 더 높았다.

마찬가지로 결론 표현의 영향이 있다, 영향을 받다는 association을 설명하는 것일 뿐, 인과관계의 의미를 내포하지 않는다.

1.5.5 시각화

(betas.year <- coef(lm(logprice~year+dong+year*dong, data = ht)))

##   (Intercept)          year        dong리   year:dong리 
## -146.35674569    0.08028525   38.07241574   -0.01917359

slope.year.d <- betas.year[2] + betas.year[4]
slope.year.l <- betas.year[2]
intercept.year.d <- betas.year[1] + betas.year[3]
intercept.year.l <- betas.year[1]

par(mfrow = c(1,3))
plot(logprice ~ year, data=ht,
            pch=20,
            xlab="year",
            ylab="log2Price", 
            main="House Price by year\n(Interaction Model)")

curve(intercept.year.d + slope.year.d*x, col="aquamarine3", add=TRUE)
curve(intercept.year.l + slope.year.l*x, col="aquamarine4", add=TRUE)

legend("topleft", c("동", "리"), 
                  lty=1, lwd=2, 
                  col=c("aquamarine3", "aquamarine4"),
                  text.col=c("aquamarine3", "aquamarine4"))

betas.area <- coef(lm(logprice~area+dong+area*dong, data = ht))

slope.area.d <- betas.area[2] + betas.area[4]
slope.area.l <- betas.area[2]
intercept.area.d <- betas.area[1] + betas.area[3]
intercept.area.l   <- betas.area[1]

plot(logprice ~ area, data=ht,
            pch=20,
            xlab="area",
            ylab="log2Price", 
            main="House Price by area\n(Interaction Model)")

curve(intercept.area.d + slope.area.d*x, col="aquamarine3", add=TRUE)
curve(intercept.area.l + slope.area.l*x, col="aquamarine4", add=TRUE)

legend("topleft", c("동", "리"), 
                  lty=1, lwd=2, 
                  col=c("aquamarine3", "aquamarine4"),
                  text.col=c("aquamarine3", "aquamarine4"))

betas.floor <- coef(lm(logprice~floor+dong+floor*dong, data = ht))

slope.floor.d <- betas.floor[2] + betas.floor[4]
slope.floor.l <- betas.floor[2]
intercept.floor.d <- betas.floor[1] + betas.floor[3]
intercept.floor.l   <- betas.floor[1]

plot(logprice ~ floor, data=ht,
            pch=20,
            xlab="floor",
            ylab="log2Price", 
            main="House Price by floor\n(Interaction Model)")

curve(intercept.floor.d + slope.floor.d*x, col="aquamarine3", add=TRUE)
curve(intercept.floor.l+ slope.floor.l*x, col="aquamarine4", add=TRUE)

legend("topleft", c("동", "리"), 
                  lty=1, lwd=2,
                  col=c("aquamarine3", "aquamarine4"),
                  text.col=c("aquamarine3", "aquamarine4"))

각각 동, 리일 때 year, area, floor의 회귀 계수 값은 달라진다. 그 중에서 특히 차이가 큰 year, logPrice 의 interaction 항에 대해 다루었고, 그래프를 봤을 때 area, floor 변수에서도 interaction 항을 찾을 여지는 보인다. Main effect model의 시각화에서와 달리동, 리 직선이 평행하게 나타난다.

1.6.1 회귀분석 가정 검토

lm_1 <- lm(price ~ year + area + floor + dong_num, data = ht)
lm_main <- lm(log2(price) ~ year + area + floor + dong_num, data = ht)
par(mfrow = c(1,2))
plot(lm_1, 1, caption = "Y = price"); plot(lm_main, 1, caption = "Y = log2price")

ht$logprice <- log2(price)

vif(lm_main)

##     year     area    floor dong_num 
## 1.071106 1.019534 1.038374 1.038475

par(mfrow = c(1,2))
plot(rstandard(lm_main), main = "log2price~."); plot(lm_main, 3)

lm_2 <- lm(log2(price) ~ year + log2(area) + floor + dong_num, data = ht)
par(mfrow = c(1,2))
plot(rstandard(lm_2), main = "log2price~log2area + ."); plot(lm_2, 3)

ncvTest(lm_main)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 10.01159, Df = 1, p = 0.0015556

price_new <- caret::BoxCoxTrans(price)
ht <- cbind(ht, price_new = predict(price_new, price))
lm_bc <- lm(price_new ~ year + area + floor + dong, data = ht)
ncvTest(lm_bc)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.5046364, Df = 1, p = 0.47747

plot(lm_main, 2)

shapiro.test(lm_main$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  lm_main$residuals
## W = 0.9874, p-value = 0.007233

lm_3 <- lm(logprice ~ year + sqrt(area) + floor + dong_num, data = ht)
plot(lm_3, 2)

shapiro.test(lm_3$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  lm_3$residuals
## W = 0.9852, p-value = 0.00236

1.6.2 가정검토 소결

종속변수에 log2 scale을 적용한 회귀 모델은 선형성과 독립성을 만족하며 등분산성과 정규성을 만족시키지 못한다. 등분산성은 box-cox transformation을 통해 해결할 수 있다.

1.6.3 Influential Point: Outlying Point, Leverage Point

설명변수 계수\(\beta_1\)의 hat value, Cook’s D, DFBETAS, RSTANDARD & RSTUDENT를 통해 확인한다.

lm_main_1 <- lm(log2(price) ~ year + area + floor + dong, data = ht_1)
ols_plot_resid_lev(lm_main_1)

ols_plot_resid_stud_fit(lm_main_1)

원 데이터에서 Outlier and Leverage point로 11, 19, 199, 215 등 4개의 점이 확인된다. 이 중에서도 Deleted Studentized Residual을 확인할 때 11, 19의 영향력이 특히 더 크고, 199, 215는 그와 비슷한 수준의 Deleted 스튜던트 잔차를 가지는 점들이 있다.

조금더 자세히 알아볼 수 있다.

influence.test <- data.frame(ht_1, cooks = cooks.distance(lm_main_1), hat = hatvalues(lm_main_1), standardized = rstandard(lm_main_1), studentized = rstudent(lm_main_1), diffinbeta = dfbeta(lm_main_1)[,2:3])

Cook’s Distance까지 포함된 influence.test라는 데이터 프레임을 생성해 Outlier, Leverage 정도가 종합적으로 큰 Influential Point를 구한다.

head(influence.test[order(influence.test$cooks, decreasing = T),])

##     dong     area price floor year dong_num logprice      cooks        hat
## 19    동  15.7838  6800    11 2013        1 12.73132 0.30792675 0.03635433
## 11    동  19.5900  6800     7 2012        1 12.73132 0.26397066 0.03357270
## 215   동 170.4840 24000    13 2001        1 14.55075 0.10939394 0.04942638
## 309   리 197.9000 27500     7 2000        0 14.74714 0.05334149 0.08738695
## 199   리  59.5900  5000     5 1991        0 12.28771 0.05157997 0.04020274
## 175   리  84.8600 22000    17 2015        0 14.42522 0.02875581 0.03577481
##     standardized studentized diffinbeta.year diffinbeta.area
## 19     -6.388360   -6.836428   -0.0012915566    6.652110e-04
## 11     -6.163882   -6.562572   -0.0011722759    6.021558e-04
## 215    -3.243360   -3.293670    0.0006639196   -4.091921e-04
## 309    -1.668928   -1.673693    0.0002691270   -2.690491e-04
## 199    -2.481347   -2.501978    0.0007451075    9.151530e-05
## 175    -1.968559   -1.977634   -0.0004978343    2.963973e-05

내림차순했을 때 Cook’s Distance가 가장 높은 두 점은 아까 Outlier & Leverage Point였던 11, 19번 점으로, 두 아파트 모두 이지역 통상적인 아파트와 비교해 면적이 현저하게 작다는 특징이 있다.(약 4~6평) 아파트명을 검색해본 결과 ‘금촌동 보광그랑베르’ ‘금촌동 그린’ 검색해본 결과 아파트가 아니라 도시형 생활주택이었으며, 이러한 작은 평수가 있고, 아파트 매매가 거래에 포함된 것으로 추정된다.

head(influence.test[order(influence.test$area, decreasing = F),])

##     dong    area price floor year dong_num logprice        cooks        hat
## 19    동 15.7838  6800    11 2013        1 12.73132 0.3079267487 0.03635433
## 11    동 19.5900  6800     7 2012        1 12.73132 0.2639706551 0.03357270
## 43    동 39.9900 10000     6 1997        1 13.28771 0.0034849491 0.02258027
## 45    동 39.9900 11000    12 1997        1 13.42522 0.0012378397 0.02290468
## 202   리 43.9200  7700    11 2000        0 12.91064 0.0244545973 0.03378433
## 44    동 49.9700 11800    10 1997        1 13.52650 0.0009656461 0.01806169
##     standardized studentized diffinbeta.year diffinbeta.area
## 19    -6.3883596  -6.8364284   -0.0012915566    6.652110e-04
## 11    -6.1638819  -6.5625723   -0.0011722759    6.021558e-04
## 43    -0.8684787  -0.8681391    0.0001743077    5.151057e-05
## 45    -0.5138348  -0.5132337    0.0001155029    3.066562e-05
## 202   -1.8700131  -1.8774930    0.0001567520    1.212534e-04
## 44    -0.5123384  -0.5117377    0.0001137176    2.315097e-05

다음으로는 Hat value를 통해 Leverage Point를 확인했다. 기울기의 분산이나 \(R^2\) 등에 영향을 미친다.

head(influence.test[order(influence.test$hat, decreasing = T),])

##     dong    area price floor year dong_num logprice       cooks        hat
## 310   리 197.900 32000    11 2000        0 14.96578 0.015828607 0.08772620
## 309   리 197.900 27500     7 2000        0 14.74714 0.053341493 0.08738695
## 215   동 170.484 24000    13 2001        1 14.55075 0.109393939 0.04942638
## 149   동 158.158 32000     1 2001        1 14.96578 0.003675971 0.04420528
## 199   리  59.590  5000     5 1991        0 12.28771 0.051579969 0.04020274
## 234   동 155.166 62500    21 2009        1 15.93157 0.002679625 0.03884859
##     standardized studentized diffinbeta.year diffinbeta.area
## 310   -0.9072026  -0.9069470    1.613725e-04   -1.460534e-04
## 309   -1.6689281  -1.6736930    2.691270e-04   -2.690491e-04
## 215   -3.2433596  -3.2936702    6.639196e-04   -4.091921e-04
## 149   -0.6304001  -0.6297961    9.351836e-05   -6.854610e-05
## 199   -2.4813469  -2.5019780    7.451075e-04    9.151530e-05
## 234    0.5757452   0.5751333   -1.238918e-05    5.658245e-05

위의 차트와 일치하는 결과를 얻을 수 있다.

1.6.4 소결

Rstandard, Rstudentized 그리고 Hat Value 등 다양한 방법을 통해 동일한 Outlier Point와 Leverage Point를 확인하고 Cook’s Distance를 통해 influential Point를 알아봤고, 공통적으로 11, 19번째 샘플을 결과에서 찾을 수 있다. 다른 아파트 매매가와 같은 선상에서 비교하기 어려운 샘플임 역시 확인했다.

1.6.5 Interval estimation

knitr::kable(confint(lm_main, level = 0.95))

2.3.1 Correlation

data01 <- data %>% dplyr::select(출석률, 결석률, 선수, 나이, 득표율)

cormat <- cor(data01, use="pairwise.complete.obs")

corrplot(cormat, diag = FALSE, method = "ellipse", type = "lower", addCoef.col="black")

ggpairs(data01,
      panel =function(x, y) {
         points (x, y)
         abline(lm(y ~ x), lty= "dashed")
          },
      diag.panel=function(x) {
        par(new=TRUE)
        hist(x, main="", axes=FALSE) 
      })

ggplot() + geom_point(data = data, aes(x=득표율, y=출석률), color='grey') + geom_smooth(data = data, aes(x=득표율, y=출석률), method = "lm", se = FALSE, color = 'black')

data001 <- data %>% filter(선수 < 3)
data002 <- data %>% filter(선수 > 3)
ggplot() + geom_point(data = data001, aes(x=득표율, y=출석률), color='darkGreen') + geom_smooth(data = data001, aes(x=득표율, y=출석률), method = "lm", se = FALSE, color = 'darkgreen') + geom_point(data = data002, aes(x=득표율, y=출석률), color='red') + geom_smooth(data = data002, aes(x=득표율, y=출석률), method = "lm", se = FALSE, color = 'red')

summary(car::powerTransform(data$출석률))

## bcPower Transformation to Normality 
##             Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## data$출석률    4.2976         4.3       3.0822       5.5131
## 
## Likelihood ratio test that transformation parameter is equal to 0
##  (log transformation)
##                            LRT df       pval
## LR test, lambda = (0) 68.08879  1 < 2.22e-16
## 
## Likelihood ratio test that no transformation is needed
##                            LRT df       pval
## LR test, lambda = (1) 36.38922  1 1.6159e-09

• 출석률(y1)은 오른쪽으로 치우쳐있고, 결석률(y2)는 왼쪽으로 치우친 형태
• 선수(x1)은 왼쪽으로 치우친 형태, 나이(x2)와 득표율(x3)은 정규분포와 유사해보임.
• 선수와 출석률은 음(-)의 상관관계, 선수와 결석률은 양(-)의 상관관계를 보인다.
• 나이도 선수와 유사하게 결석률, 출석률과 상관관계를 보이며, 나이와 선수는 높은 상관관계가 있다. 추후, 다중공선성 문제가 발생할 수 있다.
• 득표율과 출석률은 음(-)의 상관관계, 득표율과 결석률은 양(+)의 상관관계를 보인다.
• 선수는 득표율에 따른 출석률에 영향을 미치는 것으로 보인다. 따라서 선수와 득표율의 interaction term을 살펴볼 필요가 있다.

2.4.1 Multicollinearity

library(knitr)
kable(car::vif(model3))

• 초선재선항과 초선재선*득표율 항의 vif가 높은 것은 interaction term이기 때문이며, 당적 더미 변수의 vif가 높은 것은 더미 변수이기 때문이다.

2.4.2 linearity assumption

model3.1 <-  lm(출석률 ~ 득표율 + 득표율*초선재선 + 통합 + 민주 + 민생 + 정의 + 친박자유, data = dataformodel3)

par(mfrow=c(2,1))
plot(model3.1, which = 1)
plot(model3, which = 1)

• 선형성 가정을 위배하지 않는다.

2.4.3 normal distribution assumption

plot(model3.1, which = 2)

plot(model3, which = 2)

shapiro.test(resid(model3))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model3)
## W = 0.97978, p-value = 0.00194

• Q-Q plot을 보면 조금씩 잔차의 분포가 정규분포에서 떨어진 모습을 보인다. 즉 정규성 가정을 위배한다.
  
• y변수를 네제곱해줬음에도, shapiro-wilk test를 통과하지 못했다.
• 이점이 가장 어려운 문제이다.

2.4.4 Homoscedasticity Assumption

plot(model3.1, which = 3)

plot(model3, which = 3)

• 등분산성 가정을 만족한다.

bptest(model3)

## 
##  studentized Breusch-Pagan test
## 
## data:  model3
## BP = 13.562, df = 8, p-value = 0.09393

• Breusch-Pagan test에서 귀무가설이 기각되지 않아, 등분산 가정을 만족하는 것으로 한다.

2.4.5 Outlier & leverage

cutoff <- 4/(nrow(data)-length(model3$coefficients)-2)
plot(model3, which = 5, cook.levels=cutoff)
abline(h=cutoff, lty=2, col="red")

plot(model3, which = 4)

par(mfrow=c(1,1))

• cook’s distance가 큰 96번, 244번, 246번은 각각 민주당 김부겸 의원, 국회의장을 맡은 민주당 출신 문희상 의원, 국민의당 출신 무소속 김경진 의원이다.

hatleverage <- data.frame(hat = hatvalues(model3), standardized = rstandard(model3), studentized = rstudent(model3), cooks = cooks.distance(model3))


kable(head(hatleverage[order(hatleverage$cooks, decreasing = T), ]))

kable(head(hatleverage[order(hatleverage$studentized, decreasing = T), ]))

kable(head(hatleverage[order(hatleverage$standardized, decreasing = T), ]))

kable(head(hatleverage[order(hatleverage$hat, decreasing = T), ]))

• cook’s distance: 문희상 > 김경진 > 김부겸

• studentized residual: 문희상 > 정성호 > 이종구

• standardized residual: 문희상 > 정성호 > 이종구

• hat: 심상정 > 김종훈 > 조원진

• 문희상 의원은 국회의장이기 때문에 출석률이 매우 높기 때문에 outlier로 나타난다.

dataformodel4 <- dataformodel3 %>% filter(의원 != '문희상')
model4 <- lm(출석률^4 ~ 득표율 + 득표율*초선재선 + 통합 + 민주 + 민생 + 정의 + 친박자유, data = dataformodel4)
plot(model4, which = 2)

shapiro.test(resid(model4))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model4)
## W = 0.97168, p-value = 0.000126

• 문희상 의원을 제거해도 여전히 정규성을 만족하지 못한다.

outlierTest(model3)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 7 -3.511555         0.00053819      0.12648

• 본페로니 outlier 테스트에서 outlier는 나타나지 않았다.

2.4.6 interval estimation

kable(confint(model3, level = 0.95))